*****************************************************************
************* 		CREATE_TABLES_AND_FIGURES.DO 	*************
/*****************************************************************
This file loads in files created by main.m and constructs 
all of the paper's main tables and figures.

Date: January 2022
*****************************************************************/

// set paths
global main "/Users/nicolegorton/Dropbox/OTN_LAC/"
cd "${main}/generated/"

// determine folder structure 
do "${main}/code/folders"

// load in the summary stats 
import excel using "${raw_data}sumstats_Cross-Country_raw.xls",firstrow clear
replace IncomeperCapitathousandsUSD=IncomeperCapitathousandsUSD*1000

tempfile sumstats
save `sumstats'

// load in counterfactual resuts 
import excel using "${raw_data}all_cfactualsLAC.xls", firstrow clear

// start with baseline results  
keep if robust == 1
gen cross = inlist(ICC, "MERCOSUR", "ANDES")

// merge in summary statistics files
//merge m:1 Country using `roads', gen(r)
merge m:1 ICC using `sumstats', nogen 

// for each country and counterfactual, tag a single observation
bys Country cfac_type: gen first=_n==1

//define mobility as 0/1
g temp=Mobility
drop Mobility
g Mobility=1 if temp=="on"
replace Mobility=0 if temp=="off"

* order and sort
order Country ICC location cfactual_id paramset cong Mobil cfac_type
sort Country ICC location cfactual_id paramset cong Mobil cfac_type

* define some outcome variables
g y_calib = Y_calib/L_calib
g y_exp = Y_exp/L_exp

g c_calib = C_calib/L_calib
g c_exp = C_exp/L_exp

// save some unnormalized versions
gen c_calib_old = c_calib 
gen c_exp_old = c_exp

gen C_calib_old = C_calib
gen C_exp_old = C_exp

g XMY_exp =  XY_exp+ MY_exp
g XMY_calib =  XY_calib+ MY_calib

rename Population L_country
rename IncomeperCapitathousandsUSD y_country
rename AverageInfrastructure avI_country
rename TotalKMindiscretizedNetwork totKM_country


* normalize calibrated Z, H and c
g h_calib = H_calib/L_calib
	
bysort ICC Mobility gamma cfac_type: egen avZ = mean(Z_calib)
by ICC Mobility gamma cfac_type: egen avh = mean(h_calib)
replace Z_calib=Z_calib/avZ
replace h_calib=h_calib/avh
	
bysort ICC Mobility gamma cfac_type: egen avc_calib = mean(c_calib)
replace c_calib=c_calib/avc_calib
	
bysort ICC Mobility gamma cfac_type: egen avc_exp = mean(c_exp)
replace c_exp=c_exp/avc_exp

bysort ICC Mobility gamma cfac_type: egen avC_calib = mean(C_calib)
replace C_calib=C_calib/avC_calib
	
local vars = "Z_calib L_calib L_exp L_data Y_calib Y_exp Y_data C_calib c_calib C_exp y_calib y_exp I_obs I_exp c_exp XY_calib XY_exp MY_calib MY_exp XMY_exp XMY_calib"
foreach var in `vars'{
	g l`var'=log(`var')
}

local vars = "L C y c XY MY XMY Y"
foreach var in `vars'{
	g g`var' = ( `var'_exp-`var'_calib )/( 0.5*( `var'_exp+`var'_calib ) )
}

local vars =  "L C y c XY MY XMY Y"
foreach var in `vars'{
	g g`var'2 = l`var'_exp-l`var'_calib 
}

g dI = I_exp-I_obs
g gI = ( I_exp-I_obs )/( 0.5*( I_exp+I_obs ) )

* define parameter space
g low_beta = (beta==0.13)
g convex = (beta>gamma)

* keep low beta
keep if low_beta==1

* check unique id
//isid low_beta convex Mobility cfac_type ICC location
sort low_beta convex Mobility cfac_type ICC location

* tag unique cfactual
egen unique_cfac = group( low_beta convex Mobility cfac_type ICC )
egen tag_cfac = tag(unique_cfac)

egen ICCid = group(ICC)
tabulate ICC, gen(icc_)
gen logrugged = log(rugged)

bys cfactual_id ICC: egen lavR = mean(rugged)  // average c_exp within country
replace lavR = log(lavR)

label variable lc_calib "Consumption per capita"
label variable ly_calib "Income per capita"
label variable lL_calib "Population"
label variable lI_obs  "Infrastructure"
label variable diff_producer "Differentiated producer"


******************************************************************************
************* Infrasructure growth, Overall Regressions
******************************************************************************
preserve
keep if Mobility == 0
keep if cross==0

eststo clear
eststo: reghdfe gI lI_obs if cfac_type=="exp_eng", absorb(ICC) cluster(ICC)
estadd local icc "X"

eststo: reghdfe gI lI_obs  lL_calib if cfac_type=="exp_eng", absorb(ICC) cluster(ICC)
estadd local icc "X"

eststo: reghdfe gI lI_obs lL_calib ly_calib diff_ if cfac_type=="exp_eng", absorb(ICC) cluster(ICC)
estadd local icc "X"

eststo: reghdfe gI lI_obs if cfac_type=="mis_eng", absorb(ICC) cluster(ICC)
estadd local icc "X"

eststo: reghdfe gI lI_obs  lL_calib if cfac_type=="mis_eng", absorb(ICC) cluster(ICC)
estadd local icc "X"

eststo: reghdfe gI lI_obs lL_calib ly_calib diff_  if cfac_type=="mis_eng", absorb(ICC) cluster(ICC)
estadd local icc "X"

esttab using "${save_tables}regression_pooled.rtf", replace label nogap onecell ///
compress noconstant  title({\b Table 1.} {\i Outcome variable: Change in Infrastructure}) se /// r2 se ///
s(icc N r2, label("Country FE" "N" "R2")) mtitles("Expansion" "Expansion" "Expansion" ///
 "Reallocation"  "Reallocation"  "Reallocation")  
restore 

******************************************************************************
************* Figure 3: Model Fits from Calibration 
******************************************************************************
preserve

* keep benchmark
keep if Mobility == 0
keep if Ngoods==10 & city_allocation == "largest cities"
keep if convex==1
keep if cfac_type=="exp_eng"
keep if cross==0
				
sort Mobility ICC
by Mobility ICC: egen totY_calib=total(Y_calib)
by Mobility ICC: egen totY_data=total(Y_data)

replace Y_calib = Y_calib/totY_calib
replace Y_data = Y_data/totY_data

* income: model vs. data
local x = "Y_data"
local y = "Y_calib"

reg `y' `x' ,  r
local b_fixed = round(_b[`x']*1000)/1000
local se_fixed = round( _se[`x']*1000)/1000

twoway   (scatter `y' `x' , msize(vsmall) mcolor(red) ) (lfit `x' `x' , lwidth(medium) color(black) ), ///
		xtitle("Income Shares (Data)") ytitle("Income Shares (Model)") graphregion( color(white) ) ///
		legend( order( 1 2 ) lab( 1 "Fixed Labor" ) lab( 2 "45 Degree Line" ) rows(1) cols(2) pos(6)) ///
		note( "Linear regression slope (robust SE): Fixed Labor: `b_fixed' (`se_fixed')"  ) 
graph export "${images}/fit_Y.pdf", replace

restore


******************************************************************************
************* Welfare gains across countries (bar graph)
******************************************************************************
preserve
keep if Mobility == 0
keep if cross==0

keep Country cfac_type welfare_gain
duplicates drop 

replace welfare_gain = (welfare_gain-1)*100
label variable welfare_gain "Welfare gain, as % of consumption"

gen cfac_type2 = "Reallocation" if cfac_type == "mis_eng"
replace cfac_type2 = "Expansion" if cfac_type == "exp_eng"

separate welfare_gain, by(cfac_type2)
collapse (firstnm) welfare_gain1 welfare_gain2, by(Country)
graph hbar (asis) welfare_gain1 welfare_gain2, over(Country, sort(1) desc)  ///
blabel(bar, position(outside) format(%2.1f)) bar(1, bcolor(green)) bar(2, bcolor(orange)) ///
legend(pos(6) order(1 "50% Expansion" 2 "Reallocation") cols(2)) bargap(100) b1title(Welfare gain as % of consumption) ytitle("")
gr export "${images}/welfare_both.png", replace 

graph hbar (asis) welfare_gain1 , over(Country, sort(1) desc)  ///
blabel(bar, position(outside) format(%2.1f)) bar(1, bcolor(green)) bar(2, bcolor(orange)) ///
legend(off) bargap(100) b1title(Welfare gain as % of consumption) ytitle("")
gr export "${images}/welfare_exp.png", replace 

graph hbar (asis) welfare_gain2 , over(Country, sort(1) desc)  ///
blabel(bar, position(outside) format(%2.1f)) bar(1, bcolor(green)) bar(2, bcolor(orange)) ///
legend(off) bargap(100) b1title(Welfare gain as % of consumption) ytitle("")
gr export "${images}/welfare_mis.png", replace 

restore

******************************************************************************
************* Welfare gains across countries (map graph)
******************************************************************************
preserve
keep if Mobility == 0
keep if cross==0
keep Country ICC cfac_type welfare_gain L_country
duplicates drop 

replace welfare_gain = (welfare_gain-1)*100
label variable welfare_gain "Welfare gain, as % of consumption"

tabstat welfare_gain ,by(cfac_type)
tabstat welfare_gain [aw=L_country] ,by(cfac_type)

gen ISO2 = ICC
// merge on shapefile of the world 

merge m:1  ISO2 using "${shapefile}/usdb.dta", keep(match master)

spmap welfare_gain using "${shapefile}/uscoord" if cfac_type=="exp_eng", id(id) fcolor(RdYlGn) ///
osize(vvthin ...) ndsize(vvthin ...) clmethod(custom) clbreaks(0 0.25 0.5 0.75 1 1.25 1.5 2 2.5) ///
legstyle(2) legjunction(" to ")  legtitle("Welfare gain, as % of consumption")  
gr export "${images}/map_welfare_expansion.png", replace 

spmap welfare_gain using "${shapefile}/uscoord" if cfac_type=="mis_eng", id(id) fcolor(RdYlGn) ///
osize(vvthin ...) ndsize(vvthin ...) clmethod(custom) clbreaks(0 0.25 0.5 0.75 1 1.25 1.5 2 2.5) ///
legstyle(2) legjunction(" to ")  legtitle("Welfare gain, as % of consumption") 
gr export "${images}/map_welfare_reallocation.png", replace 
restore 

******************************************************************************
************* Comparison to World Bank Projects 
******************************************************************************
preserve
import delimited using "${raw_data}wb_Iexp.csv", clear
drop if var1 == "Venezuela" | var1 == "Panama"

graph hbar (asis) corrs_amount , over(var1, sort(1) desc)  ///
blabel(bar, position(outside) format(%3.2f)) bar(1, bcolor(blue)) bar(2, bcolor(red)) ///
legend(off)   exclude0 ///
 bargap(100) b1title("Correlation Between Infrastructure Growth under Expansion and World Bank Projects ($)", size(small)) ytitle("") ylab(-0.2(0.2)1) 
gr export "${images}/wb_amount.png", replace 

 graph hbar (asis) corrs_count , over(var1, sort(1) desc)  ///
blabel(bar, position(outside) format(%3.2f)) bar(1, bcolor(blue)) bar(2, bcolor(red)) ///
legend(off)  exclude0   ///
bargap(100) b1title("Correlation Between Infrastructure Growth under Expansion and World Bank Projects (#)", size(small)) ytitle("") ylab(-0.2(0.2)1) 
gr export "${images}/wb_count.png", replace 
restore 

******************************************************************************
************* Maps and scatterplots 
******************************************************************************
** Scatterplots 
preserve 
keep if Mobil==0 
drop if cross==1

drop gC 
g gC = ( C_exp_old-C_calib_old )/( 0.5*(C_exp_old+C_calib_old) ) * 100

// all countries together 
local x = "gC"
local y = "ly_calib"
reg `y' `x' if cfac_type=="exp_eng", r
local b = round(_b[`x']*1000)/1000
local se = round( _se[`x']*1000)/1000
			
graph twoway (scatter gC ly_calib if  cfac_type=="exp_eng" & !inlist(Country, "Brazil", "Argentina", "Mexico"), mcolor(black) msymbol(Oh)) ///
(scatter gC ly_calib if  cfac_type=="exp_eng" & inlist(Country, "Brazil"), mcolor(blue) msymbol(D)) ///
(scatter gC ly_calib if  cfac_type=="exp_eng" & inlist(Country, "Argentina"), mcolor(red) msymbol(o)) ///
(scatter gC ly_calib if  cfac_type=="exp_eng" & inlist(Country, "Mexico"), mcolor(green) msymbol(X)) ///
(lfit gC ly_calib  if cfac_type=="exp_eng", lcolor(blue) lpattern(solid)) , ///
note( "Linear regression slope (robust SE): `b' (`se')", size(small)) ///
xtitle("log(Initial Income per Capita)") ytitle("Percentage Change in Consumption") ///
legend( order(4 "Mexico" 3 "Argentina" 2 "Brazil" 1 "All Other Countries") cols(4))
gr export "${images}/all_correlation.png", replace

restore 

** Maps, generated by 
preserve 
cd "${rawdata}gridmaps"
global countries "${raw_data_maps}COUNTRY_SHP/countries_shp/"

local countries Argentina Brazil

foreach c in `countries' {
	cap shp2dta using `c', database(`c'_db) coordinates(`c'_coords) genid(id)
	cap shp2dta using `c'_roads, database(`c'_roads_db) coordinates(`c'_roads_coords) genid(id)
	
	use `c'_roads_coords, clear
	gen id = _ID
	merge m:1 id using `c'_roads_db, nogen
	keep if dI > 0.001 & dI<.
	
	gen percentile = .
	sum dI , det
	replace percentile = 0 if dI<r(p50)
	replace percentile = 50 if dI>=r(p50) &  dI<r(p75)
	replace percentile = 75 if dI>=r(p75) &  dI<r(p90)
	replace percentile = 90 if dI>=r(p90) 
			
	save `c'_roads_coords_wI,replace
	
	use `c'_db, clear

	sum chat if chat<1, det
	local minlow=r(min)
	local p25low =r(p25) 
	local p75low =r(p75) 
	local maxlow=r(max)

	sum chat if chat>1, det
	local minhi=r(min)
	local p25hi =r(p25) 
	local p75hi =r(p75) 
	local maxhi= r(max)

	di `maxhi'
	di `minlow' `maxlow' 0.9999999 1 1.0000001 `minhi' `maxhi'

	if "`c'" != "Peru" {
		spmap chat using `c'_coords, id(id) ///
		fcolor(orange orange*.7 orange*.4 orange*.1  white green*.1 green*.4 green*.7 green) ocolor(none ...) ///
		clbreaks(`minlow' `p25low' `p75low' `maxlow' 0.9999999 1 `minhi' `p25hi' `p75hi' `maxhi') clmethod(custom) ///
		legjunction("to") legorder(lohi) legend(label(2 "{bf:Decrease}") label(3 " ") label(4 " ") label(5 " ") ///
		label(6 "{bf:No change}" ) label(7 " ") label(8 " ") label(9 " ")  label(10 "{bf:Increase}") row(1) pos(6))   ///
		legtitle("{bf:Change:} ") ///
		line(data("`c'_roads_coords_wI")  ///
		color(black*.3 black*.4 black*.7 black) size(0.2 0.4 0.6 1.5 ) legenda(off)  by(percentile))
	}
	
	if "`c'" == "Peru" {
		
		spmap chat using `c'_coords, id(id) ///
		fcolor(orange orange*.7 orange*.4 orange*.1  white green*.1 green*.4 green*.7 green) ocolor(none ...) ///
		clbreaks(`minlow' `p25low' `p75low' `maxlow' 0.9999999 1 `minhi' `p25hi' `p75hi' `maxhi') clmethod(custom) ///
		legjunction("to") legorder(lohi) legend(label(2 "{bf:Decrease}") label(3 " ") label(4 " ") label(5 " ") ///
		label(6 "{bf:No change}" ) label(7 " ") label(8 " ") label(9 " ")  label(10 "{bf:Increase}") row(1) pos(6))   ///
		legtitle("{bf:Change:} ") ///
		line(data("`c'_roads_coords_wI")  ///
		color(black*.3 black*.4 black*.7 black) size(0.2 0.4 0.6 1.5 ) legenda(off)  by(percentile)) ///
		polygon(data(${countries}/uscoord) select(keep if _ID==188))
		
		
	}
	
	gr export `c'.png, replace
}
restore 


******************************************************************************
************* Labor mobility robustness check (correlation)
******************************************************************************
preserve
keep if cross==0
keep ICC Mobil welfare_gain cfac_type
drop if ICC=="VE" // the model never calibrated and the fit is terrible for income (like, really bad)
duplicates drop 
replace welfare_gain = (welfare_gain-1)*100
reshape wide welfare_gain, i(ICC cfac_) j(Mobility)

correlate welfare_gain1 welfare_gain0
local rho = round(`r(rho)'*100, 1) 

di `rho'

graph twoway (scatter welfare_gain1 welfare_gain0 if cfac=="exp_eng", mlabel(ICC) msymbol(0) mcolor(red) mlabcolor(red)) ///
 (scatter welfare_gain1 welfare_gain0 if cfac=="mis_eng", mlabel(ICC) msymbol(o) mcolor(blue) mlabcolor(blue))  ///
 (line welfare_gain0 welfare_gain0 if cfac=="mis_eng", lpattern(solid)), ///
 legend(order(1 "Expansion" 2 "Reallocation" 3 "45 Degree Line") pos(6) cols(3)) ///
 xtitle("Welfare Gain, Fixed Labor") ytitle("Welfare Gain, Mobile Labor") ///
 note("Correlation: `rho'%")
gr export "${images}welfare_mobil_correlation.png", replace
restore 

******************************************************************************
************* Country-level Summary Statistics
******************************************************************************
** exported directly from Matlab and Python files, see readme

******************************************************************************
************* Labor mobility, regression table 
******************************************************************************
// regression to understand population mobility 
preserve 
drop if ICC=="VE" // the model never calibrated and the fit is terrible for income (like, really bad)
keep if Mobility == 1

label variable lc_calib "Consumption per capita"
label variable ly_calib "Tradable Income per capita"
label variable lL_calib "Population"
label variable lI_obs  "Infrastructure"
label variable gI "Infrastructure Growth"
label variable diff_producer "Differentiated Producer"

eststo clear
eststo: reg gI lI_obs lL_calib ly_calib diff if cfac_type=="exp_eng",  absorb(ICC) cluster(ICC)
eststo: reg gI lI_obs lL_calib ly_calib diff if cfac_type=="mis_eng",  absorb(ICC) cluster(ICC)
eststo: reg gL lL_calib ly_calib lc_calib lI_obs gI diff_producer if cfac_type=="exp_eng",  absorb(ICC) cluster(ICC)
eststo: reg gL lL_calib ly_calib lc_calib lI_obs gI diff_producer if cfac_type=="mis_eng",  absorb(ICC) cluster(ICC)

esttab using "${save_tables}regression_pooled_mobility.rtf", replace label nogap onecell ///
compress noconstant  title({\b Table 3.} {\i Outcome variable: Changes in Infrastructure and Population}) se /// r2 se ///
s(icc N r2, label("Country FE" "N" "R2")) mtitles("Expansion" "Reallocation" "Expansion" "Reallocation") ///
mgroups("Infrastructure Growth" "Population Growth", pattern(1 0 1 0) span)  
restore 

******************************************************************************
************* Transnational welfare gains by country
******************************************************************************
preserve // gen welfare gains by grouping 
replace welfare_gain = (welfare_gain-1)*100

tabstat welfare_gain if cross==1 & ICC=="ANDES", by(cfac_t)
tabstat welfare_gain if cross==1 & ICC=="MERCOSUR", by(cfac_t)

restore 

preserve 
keep if cross==1

gen country = ""
replace country = "Peru" if inrange(location, 1, 24) & ICC=="ANDES"
replace country = "Colombia" if inrange(location, 25, 56 ) & ICC=="ANDES"
replace country = "Ecuador" if inrange(location, 57, 80) & ICC=="ANDES"
replace country = "Bolivia" if inrange(location , 80, 88) & ICC=="ANDES"

replace country = "Argentina" if inrange(location, 1, 23) & ICC=="MERCOSUR"
replace country = "Brazil" if inrange(location, 24, 50 ) & ICC=="MERCOSUR"
replace country = "Paraguay" if inrange(location,51, 53) & ICC=="MERCOSUR"
replace country = "Uruguay" if inrange(location, 54 ,57) & ICC=="MERCOSUR"

drop c_calib c_exp
gen c_calib = c_calib_old
gen c_exp = c_exp_old

gen alpha = 0.4 // parameter for utility 

// compute welfare gains by grid cell, and then by country
gen u_calib= (c_calib/alpha)^alpha
gen u_exp = (c_exp/alpha)^alpha

bys ICC cfac_type: egen u_total_calib = total(L_calib*u_calib)
bys ICC cfac_type: egen u_total_exp = total(L_calib*u_exp)

// this recovers the total welfare change when alpha=0.4
// veryify that it matches reported figures and it does 
gen change = u_total_exp/u_total_calib
replace change= change^(1/alpha)

// compute each country's welfare gain:
bys ICC cfac_type country: egen u_total_calib_c = total(L_calib*u_calib)
bys ICC cfac_type country: egen u_total_exp_c = total(L_calib*u_exp)

// this recovers the total welfare change when alpha=0.4
gen change_c = u_total_exp_c/u_total_calib_c
replace change_c= change_c^(1/alpha)

gen counterfactual = "Expansion" if cfac_type=="exp_eng"
replace counterfactual = "Reallocation" if cfac_type=="mis_eng"

keep change_c ICC counterfactual country 
duplicates drop
replace change_c = (change_c - 1)*100

g change_c_label = round(change_c*100)/100
label variable change_c_label "Welfare Gain as % of Consumption"

graph hbar (asis) change_c_label if ICC=="MERCOSUR", over(country) over(counterfactual) ///
bar(1, fcolor(gray) lcolor(white)) blabel(total) yscale(range(-0.4 3))
gr export "${save_figures}welfare_by_country_MERCOSUR.pdf", replace

graph hbar (asis) change_c_label if ICC=="ANDES", over(country) over(counterfactual) ///
bar(1, fcolor(gray) lcolor(white)) blabel(total) 
gr export "${save_figures}welfare_by_country_ANDES.pdf", replace
restore 

******************************************************************************
************* Other
******************************************************************************
preserve
import delimited using "${raw_data}wb_Iexp.csv", clear
rename var1 Country

drop if corrs_amount == 0
label variable corrs_amount	"Corr(WB Project Spending ($), Infrastructure Growth)"

** Figure 13
graph hbar (asis) corrs_amount corrs_count, over(Country, sort(1) descending)  bar(1, fcolor(red))  bar(2, fcolor(blue)) ///
 legend(pos(6) cols(1) ///
 order(1 "Corr(WB Project Spending ($), Infrastructure Growth)" ///
 2 "Corr(# of WB Projects, Infrastructure Growth)" ))
   //bar(3, fcolor(yellow)) 
gr export "${save_figures}wbp_correlations.png", replace

** Figure A2
graph hbar (asis) corrs_population corrs_population2, over(Country, sort(1) descending)  bar(1, fcolor(green) lcolor(none))  bar(2, fcolor(purple) lcolor(none)) ///
 legend(pos(6) cols(1) ///
 order(1 "Corr(WB Project Spending ($), Population)" ///
 2 "Corr(# of WB Projects, Population)" ))
   //bar(3, fcolor(yellow)) 
gr export "${save_figures}wbp_population_correlations.png", replace

** Figure A3
graph twoway (scatter corrs_amount corrs_population, mlabel(Country))(scatter corrs_count corrs_population2, mlabel(Country)), ///
 legend(pos(6) cols(2) ///
 order(1 "WB Project Spending ($)" ///
 2 "# of WB Projects" )) ///
xtitle("Corr(WB Projects, Population)") ytitle("Corr(WB Projects, Infrastructure Growth)")
gr export "${save_figures}wbp_population_correlations2.png", replace
restore 

******************************************************************************
************* Robustness Check Figures and Scatterplots 
******************************************************************************
// load in counterfactual resuts 
import excel using "${raw_data}all_cfactualsLAC.xls", firstrow clear
drop if inlist(ICC, "MX", "AR", "CO", "PE", "BR") // focus on smaller countries 

// start with baseline results  
gen cross = inlist(ICC, "MERCOSUR", "ANDES")
drop if cross==1
keep if Mobil=="off"

replace welfare_gain = (welfare_gain-1)*100

// this constructs the table of welfare gains across robustness checks
preserve 
keep welfare_gain Country robust cfac_t
format %2.1g welfare_gain 
duplicates drop 
reshape wide welfare_gain, i(Country cfac_t) j(robust)
reshape wide welfare_gain1 welfare_gain2 welfare_gain3 welfare_gain4, i(Country) j( cfac_t, string)
order Country welfare_gain1exp_eng welfare_gain2exp_eng  welfare_gain3exp_eng welfare_gain4exp_eng ///
welfare_gain1mis_eng welfare_gain2mis_eng  welfare_gain3mis_eng welfare_gain4mis_eng

export delimited using "${save_tables}welfare_gains_by_robust.csv", replace
restore 

// this constructs the correlation of infrastructure placement 
preserve 
g dI = I_exp-I_obs
bys ICC cfac_type robust: egen avgI = mean(dI)
replace dI=dI/avgI
keep robust dI ICC location cfac_type robust
reshape wide dI, j(robust) i(ICC location cfac_type)

corr dI2 dI1
local r2 = round(`r(rho)'*100)

corr dI3 dI1
local r3 = round(`r(rho)'*100)

corr dI4 dI1
local r4 = round(`r(rho)'*100)


graph twoway (scatter dI2 dI1, msymbol(D) mcolor(blue)) ///
(scatter dI3 dI1, msymbol(0h) mcolor(red))(scatter dI4 dI1, msymbol(o) mcolor(purple))(line dI1 dI1), ///
legend(order(1 "WorldPop (corr=`r2'%)" 2 "OpenStreetMap (corr=`r3'%)" 3 "WorldPop + OpenStreetMap (corr=`r4'%)" 4 "X-Y Line") cols(2) size(small)) ///
xtitle("Baseline Infrastructure Change") ytitle("Robustness Infrastructure Change") 
gr export "${save_figures}robustness_infrastructure_correlations.png", replace
restore 





